        
function output = fun_optimize_wEC(Z,kappa,sigma,varepsilon,beta, ...
                                      R_base, t, T, MEC, ...
                                      p, s_vector, w)
        
        %-------------------------------------------------------------%
        % Needed inputs
        %-------------------------------------------------------------%
        % preferences, scalars:     varepsilon, sigma, kappa, beta
        % scalars                   R_base, t, T, MEC,
        % vectors:                  p(1:T); s_vector(1:T); w(1:T)
        %-------------------------------------------------------------%
        R = R_base;
        c = NaN(T,1);
        g = NaN(T,1);
        %----------------------------------------%
        % Solve for unshocked gt; no entry cost
        %----------------------------------------%
        
        A = (p(t)/kappa)^sigma * sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t))  );
        B = sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        C = p(t)^varepsilon * sum( p(t:T).^(1-varepsilon) .*   beta.^(varepsilon.*(s_vector(t:T)-t)) .* R.^((varepsilon-1).*(s_vector(t:T)-t))  );
      
        syms x
        gt_base = vpasolve(x^(sigma/varepsilon)*A - B + x*C == 0,x,[0.000000001 Inf]);
        

        % Now compute vectors g and c from FOCs
        g(t:T) = (p(t)./p(t:T)).^varepsilon .* beta.^((s_vector(t:T)-t)*varepsilon) .*R.^((s_vector(t:T)-t)*varepsilon)*gt_base;
       
        c(t:T) = (p(t:T)./kappa).^sigma .* g(t:T).^(sigma/varepsilon);

        % Calculate life-time utility with giving
        U = sum( beta.^(s_vector(t:T)-t).*(  c(t:T).^(1-1/sigma)./(1-1/sigma) + kappa * (g(t:T)).^(1-1/varepsilon)./(1-1/varepsilon) ) );
        
        gt_base_noEC = gt_base;
        g_base_noEC = g;
        %-----------------------------------------------------%
        % Solve for unshocked gt; constraining g=0 forever
        % c_0 * A2 = B2
        % B2 now accounts for a bonus: you get to turn MEC into consumption
        %-----------------------------------------------------------------%
        gt_base = 0;
        A2 = sum( beta.^(sigma.*(s_vector(t:T)-t)) .* R.^((sigma-1).*(s_vector(t:T)-t)));
        B = MEC + sum( w(t:T) .* R.^(-(s_vector(t:T)-t)) );
        ct_base = B/A2;
        
        g(t:T) = 0;
        for s=t+1:T
            c(s) = (beta*R)^sigma * ct_base;
        end
        
        % Calculate life-time utility with giving
        perperiod_givingutility_loss = g_base_noEC(t:T).*kappa .* (g_base_noEC(t:T)).^(-1/varepsilon); % from 1st order taylor approx
        new_perperiod_givingutility = kappa * (g_base_noEC(t:T)).^(1-1/varepsilon)./(1-1/varepsilon) - perperiod_givingutility_loss;
        U_notgive = sum(  beta.^(s_vector(t:T)-t).*(c(t:T).^(1-1/sigma)./(1-1/sigma)) + new_perperiod_givingutility );


        % Choose to give if U > U_give
        if U > U_notgive
            gt_base_wEC  = gt_base_noEC;
        elseif U<= U_notgive 
            gt_base_wEC  = 0;
        end

        UtilDiff = U-U_notgive;

        if Z == 1
            output = gt_base_wEC;
        elseif Z == 2
            output = UtilDiff;
        elseif Z==3 
            output = g_base_noEC; % giving if no EC
        end


        if isnan(output)
            error('nan....')
        end

